Python 練功房¶

# 加入函式庫和函數
from IPython.display import Audio
from numpy import pi,sin,exp,cos
import numpy as np

# 定義音階
A,B, C, D, E, F, G, A5, B5, C5,D5,E5,F5,G5,A6,B6, C6 = 57,59,\
          60, 62, 64, 65, 67, 69, 71, 72,74,76,77,79,81,83,84
scale = [C, D, E, F, G, A, B, C5 ,D5, E5,F5,G5,A6,B6,C6 ]

def sine_env(t):
    return sin(pi*t/t[-1]) #t[-1] = last value

def flute_env(t):
    return sin(pi*t/t[-1])**0.4

# 歌曲
train=[(D5,1),(E5,1), (D5,1/2),(E5,1/2),(D5,1/2),(C5,1/2),
       (D5,1/2),(E5,1/4),(E5,1/4),(D5,1/2),(C5,1/2),(D5,1),(C5,1/2),(A5,1/2),
       (G,1/2),(G,1/2),(A5,1/2),(C5,1/2),(G,2),
       (A5,1),(D5,1),(A5,1/2),(C5,1/2),(C5,1/2),(A,1/2),
       (D5,1/2),(C5,1/4),(D5,1/2),(C5,1/4),
       (A5,1/2),(C5,1/2),(A5,1/2),(C5,1/2),
       (G,1/2),(G,1/4),(G,1/4),(G,1/2),(G,1/2),
       (A5,1/2),(C5,1/2),(A5,1/2),(G,1/2),(D,3/2),(A5,1/2),(G,2)]

# Solo by flute
BaseNote = 466.16
Fs = 44100
sound = []
for note in train:
    fnum, dur = note
    t = np.linspace(0,dur/2,dur/2*Fs)
    f = BaseNote * 2**((fnum-69)/12)
    sinusoid = sin(2*pi*f*t)
    #sound.append(guitar_env(t) * sinusoid)
    sound.append(flute_env(t) * sinusoid)
train_guitar_sound = np.concatenate(sound)
Audio(train_guitar_sound, rate=Fs)

# 加入製圖函式庫
%matplotlib inline
from scipy.integrate import ode,odeint
import matplotlib.pylab as plt
from matplotlib import animation
from JSAnimation import IPython_display

# 定義參數
F=0
R=0.2
h0=2

theta0=0
v0=5.

t = np.linspace(0, 20, 2000)

# 時間軌跡計算
def dx(x, t):
    """
    The right-hand side of the pendulum ODE
    (v,theta,x,y)
    """
    x1, x2, x3, x4 = x[0], x[1], x[2], x[3]
    
    dx1 = -sin(x2)-R*x1*x1+F
    dx2 = (x1*x1-cos(x2))/x1
    
    dx3 = x1*cos(x2)
    dx4 = x1*sin(x2)
    
    return [dx1, dx2, dx3, dx4]

# 初始值
theta0Rad=theta0/180.*pi
x0 = [v0,theta0Rad,0,h0]

# 計算飛行軌跡
x = odeint(dx, x0, t)

# 計算何時落地
distanceI=(x[:,3]>0)
y=x[:,3].clip(0)
TavelTime=sum(distanceI)

# 動畫
theta0=30
v0=5.

ratio=20
anim = plt.figure(figsize=(10,3))
ax = anim.add_subplot(111)
ax.set_title("Paper Plane",fontsize=14)

plt.ylim([-1,h0+3])
plt.xlim([-5,15])
plt.plot([0,0],[0,h0],lw=8)
plt.plot([-5,15],[0,0],'--b')
plt.title('Paper Plane', size=24)
plt.text(1.5,4,'Initial velocity %2.2f, attack angle %2.2f degree' %(v0,theta0), size=16, color='grey')
plt.text(1.5,3.4, 'Traveling Distance: %2.2f' %x[TavelTime,2], size=16, color='grey')
def init():
    return None
def animate(i):
    if (i>0):
        dx=(x[i*ratio,2]-x[(i-1)*ratio,2])#/np.abs(x[i*ratio,2]-x[(i-1)*ratio,2])*0.4
        dy=(y[i*ratio]-y[(i-1)*ratio])#/np.abs(y[i*ratio]-y[(i-1)*ratio])*0.4
    else:
        dx=np.cos(theta0)*0.6;
        dy=np.sin(theta0)*0.6;
    return ax.plot(x[0:ratio*i:ratio,2],y[0:ratio*i:ratio],'--r'), \
        ax.arrow(x[i*ratio,2],y[i*ratio],dx,dy,color='red',head_width=0.4, head_length=0.4,alpha=0.3)

animation.FuncAnimation(anim, animate, init_func=init,frames=int(TavelTime/ratio))     

# 補充額外的函式庫和函數
from numpy.linalg import norm
from matplotlib import animation
from ipywidgets import interact, FloatSlider

# 參數
Mpos=np.array([[1,0]])
Epos=np.array([[0,0]])
p=29.5
c=0.15
dt=0.1
i=1

# 月亮和太空船的軌跡計算
def MoonTraj(vec):
    c=vec
    Mpos=np.array([[1,0]])
    Epos=np.array([[0,0]])

    i=1
    #while Dpos[-1][1]<Mpos[0][1] or i<40:
    while i<800:
          Mpos1=np.array([cos(2*pi*i*dt/p),sin(2*pi*i*dt/p)])  
          Evel=c*(Mpos[-1]-Epos[-1])/norm(Mpos[-1]-Epos[-1])

          Epos2 = Epos[-1]+Evel*dt

          Mpos=np.vstack([Mpos,Mpos1])
          Epos=np.vstack([Epos,Epos2]) 
          i=i+1

    figM = plt.figure(figsize=(6,6))
    #ax = figM.add_subplot(111)
    plt.title("Moon Landing",size=14)

    plt.xlim([-1.1,1.1])
    plt.ylim([-1.1,1.1])
    plt.scatter(Mpos[:,0],Mpos[:,1], c='r' )
    plt.scatter(Epos[::4,0],Epos[::4,1],c='g')            
    plt.text(-0.5, -0.5,
            "Velocity of Spacecraft={0:.2f}".format(c),
            size=14, color='gray')       
            
    #return figM

# 改變參數觀察太空登陸艇 (綠色軌跡) 如何才能登上月球(外圍的軌跡)
interact(MoonTraj,vec=FloatSlider(min=0.05, max=0.30, step=0.01))